home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / exp.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  15.9 KB  |  611 lines
  1. /* specfunc/exp.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_sf_gamma.h>
  26. #include <gsl/gsl_sf_exp.h>
  27.  
  28. #include "error.h"
  29.  
  30. /* Evaluate the continued fraction for exprel.
  31.  * [Abramowitz+Stegun, 4.2.41]
  32.  */
  33. static
  34. int
  35. exprel_n_CF(const int N, const double x, gsl_sf_result * result)
  36. {
  37.   const double RECUR_BIG = GSL_SQRT_DBL_MAX;
  38.   const int maxiter = 5000;
  39.   int n = 1;
  40.   double Anm2 = 1.0;
  41.   double Bnm2 = 0.0;
  42.   double Anm1 = 0.0;
  43.   double Bnm1 = 1.0;
  44.   double a1 = 1.0;
  45.   double b1 = 1.0;
  46.   double a2 = -x;
  47.   double b2 = N+1;
  48.   double an, bn;
  49.  
  50.   double fn;
  51.  
  52.   double An = b1*Anm1 + a1*Anm2;   /* A1 */
  53.   double Bn = b1*Bnm1 + a1*Bnm2;   /* B1 */
  54.   
  55.   /* One explicit step, before we get to the main pattern. */
  56.   n++;
  57.   Anm2 = Anm1;
  58.   Bnm2 = Bnm1;
  59.   Anm1 = An;
  60.   Bnm1 = Bn;
  61.   An = b2*Anm1 + a2*Anm2;   /* A2 */
  62.   Bn = b2*Bnm1 + a2*Bnm2;   /* B2 */
  63.  
  64.   fn = An/Bn;
  65.  
  66.   while(n < maxiter) {
  67.     double old_fn;
  68.     double del;
  69.     n++;
  70.     Anm2 = Anm1;
  71.     Bnm2 = Bnm1;
  72.     Anm1 = An;
  73.     Bnm1 = Bn;
  74.     an = ( GSL_IS_ODD(n) ? ((n-1)/2)*x : -(N+(n/2)-1)*x );
  75.     bn = N + n - 1;
  76.     An = bn*Anm1 + an*Anm2;
  77.     Bn = bn*Bnm1 + an*Bnm2;
  78.  
  79.     if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
  80.       An /= RECUR_BIG;
  81.       Bn /= RECUR_BIG;
  82.       Anm1 /= RECUR_BIG;
  83.       Bnm1 /= RECUR_BIG;
  84.       Anm2 /= RECUR_BIG;
  85.       Bnm2 /= RECUR_BIG;
  86.     }
  87.  
  88.     old_fn = fn;
  89.     fn = An/Bn;
  90.     del = old_fn/fn;
  91.     
  92.     if(fabs(del - 1.0) < 2.0*GSL_DBL_EPSILON) break;
  93.   }
  94.  
  95.   result->val = fn;
  96.   result->err = 2.0*(n+1.0)*GSL_DBL_EPSILON*fabs(fn);
  97.  
  98.   if(n == maxiter)
  99.     GSL_ERROR ("error", GSL_EMAXITER);
  100.   else
  101.     return GSL_SUCCESS;
  102. }
  103.  
  104.  
  105. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  106.  
  107. #ifndef HIDE_INLINE_STATIC
  108. int gsl_sf_exp_e(const double x, gsl_sf_result * result)
  109. {
  110.   if(x > GSL_LOG_DBL_MAX) {
  111.     OVERFLOW_ERROR(result);
  112.   }
  113.   else if(x < GSL_LOG_DBL_MIN) {
  114.     UNDERFLOW_ERROR(result);
  115.   }
  116.   else {
  117.     result->val = exp(x);
  118.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  119.     return GSL_SUCCESS;
  120.   }
  121. }
  122. #endif
  123.  
  124. int gsl_sf_exp_e10_e(const double x, gsl_sf_result_e10 * result)
  125. {
  126.   if(x > INT_MAX-1) {
  127.     OVERFLOW_ERROR_E10(result);
  128.   }
  129.   else if(x < INT_MIN+1) {
  130.     UNDERFLOW_ERROR_E10(result);
  131.   }
  132.   else {
  133.     const int N = (int) floor(x/M_LN10);
  134.     result->val = exp(x-N*M_LN10);
  135.     result->err = 2.0 * (fabs(x)+1.0) * GSL_DBL_EPSILON * fabs(result->val);
  136.     result->e10 = N;
  137.     return GSL_SUCCESS;
  138.   }
  139. }
  140.  
  141.  
  142. int gsl_sf_exp_mult_e(const double x, const double y, gsl_sf_result * result)
  143. {
  144.   const double ay  = fabs(y);
  145.  
  146.   if(y == 0.0) {
  147.     result->val = 0.0;
  148.     result->err = 0.0;
  149.     return GSL_SUCCESS;
  150.   }
  151.   else if(   ( x < 0.5*GSL_LOG_DBL_MAX   &&   x > 0.5*GSL_LOG_DBL_MIN)
  152.           && (ay < 0.8*GSL_SQRT_DBL_MAX  &&  ay > 1.2*GSL_SQRT_DBL_MIN)
  153.     ) {
  154.     const double ex = exp(x);
  155.     result->val = y * ex;
  156.     result->err = (2.0 + fabs(x)) * GSL_DBL_EPSILON * fabs(result->val);
  157.     return GSL_SUCCESS;
  158.   }
  159.   else {
  160.     const double ly  = log(ay);
  161.     const double lnr = x + ly;
  162.  
  163.     if(lnr > GSL_LOG_DBL_MAX - 0.01) {
  164.       OVERFLOW_ERROR(result);
  165.     }
  166.     else if(lnr < GSL_LOG_DBL_MIN + 0.01) {
  167.       UNDERFLOW_ERROR(result);
  168.     }
  169.     else {
  170.       const double sy   = GSL_SIGN(y);
  171.       const double M    = floor(x);
  172.       const double N    = floor(ly);
  173.       const double a    = x  - M;
  174.       const double b    = ly - N;
  175.       const double berr = 2.0 * GSL_DBL_EPSILON * (fabs(ly) + fabs(N));
  176.       result->val  = sy * exp(M+N) * exp(a+b);
  177.       result->err  = berr * fabs(result->val);
  178.       result->err += 2.0 * GSL_DBL_EPSILON * (M + N + 1.0) * fabs(result->val);
  179.       return GSL_SUCCESS;
  180.     }
  181.   }
  182. }
  183.  
  184.  
  185. int gsl_sf_exp_mult_e10_e(const double x, const double y, gsl_sf_result_e10 * result)
  186. {
  187.   const double ay  = fabs(y);
  188.  
  189.   if(y == 0.0) {
  190.     result->val = 0.0;
  191.     result->err = 0.0;
  192.     result->e10 = 0;
  193.     return GSL_SUCCESS;
  194.   }
  195.   else if(   ( x < 0.5*GSL_LOG_DBL_MAX   &&   x > 0.5*GSL_LOG_DBL_MIN)
  196.           && (ay < 0.8*GSL_SQRT_DBL_MAX  &&  ay > 1.2*GSL_SQRT_DBL_MIN)
  197.     ) {
  198.     const double ex = exp(x);
  199.     result->val = y * ex;
  200.     result->err = (2.0 + fabs(x)) * GSL_DBL_EPSILON * fabs(result->val);
  201.     result->e10 = 0;
  202.     return GSL_SUCCESS;
  203.   }
  204.   else {
  205.     const double ly  = log(ay);
  206.     const double l10_val = (x + ly)/M_LN10;
  207.  
  208.     if(l10_val > INT_MAX-1) {
  209.       OVERFLOW_ERROR_E10(result);
  210.     }
  211.     else if(l10_val < INT_MIN+1) {
  212.       UNDERFLOW_ERROR_E10(result);
  213.     }
  214.     else {
  215.       const double sy  = GSL_SIGN(y);
  216.       const int    N   = (int) floor(l10_val);
  217.       const double arg_val = (l10_val - N) * M_LN10;
  218.       const double arg_err = 2.0 * GSL_DBL_EPSILON * fabs(ly);
  219.  
  220.       result->val  = sy * exp(arg_val);
  221.       result->err  = arg_err * fabs(result->val);
  222.       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  223.       result->e10 = N;
  224.  
  225.       return GSL_SUCCESS;
  226.     }
  227.   }
  228. }
  229.  
  230.  
  231. int gsl_sf_exp_mult_err_e(const double x, const double dx,
  232.                              const double y, const double dy,
  233.                              gsl_sf_result * result)
  234. {
  235.   const double ay  = fabs(y);
  236.  
  237.   if(y == 0.0) {
  238.     result->val = 0.0;
  239.     result->err = fabs(dy * exp(x));
  240.     return GSL_SUCCESS;
  241.   }
  242.   else if(   ( x < 0.5*GSL_LOG_DBL_MAX   &&   x > 0.5*GSL_LOG_DBL_MIN)
  243.           && (ay < 0.8*GSL_SQRT_DBL_MAX  &&  ay > 1.2*GSL_SQRT_DBL_MIN)
  244.     ) {
  245.     double ex = exp(x);
  246.     result->val  = y * ex;
  247.     result->err  = ex * (fabs(dy) + fabs(y*dx));
  248.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  249.     return GSL_SUCCESS;
  250.   }
  251.   else {
  252.     const double ly  = log(ay);
  253.     const double lnr = x + ly;
  254.  
  255.     if(lnr > GSL_LOG_DBL_MAX - 0.01) {
  256.       OVERFLOW_ERROR(result);
  257.     }
  258.     else if(lnr < GSL_LOG_DBL_MIN + 0.01) {
  259.       UNDERFLOW_ERROR(result);
  260.     }
  261.     else {
  262.       const double sy  = GSL_SIGN(y);
  263.       const double M   = floor(x);
  264.       const double N   = floor(ly);
  265.       const double a   = x  - M;
  266.       const double b   = ly - N;
  267.       const double eMN = exp(M+N);
  268.       const double eab = exp(a+b);
  269.       result->val  = sy * eMN * eab;
  270.       result->err  = eMN * eab * 2.0*GSL_DBL_EPSILON;
  271.       result->err += eMN * eab * fabs(dy/y);
  272.       result->err += eMN * eab * fabs(dx);
  273.       return GSL_SUCCESS;
  274.     }
  275.   }
  276. }
  277.  
  278.  
  279. int gsl_sf_exp_mult_err_e10_e(const double x, const double dx,
  280.                              const double y, const double dy,
  281.                              gsl_sf_result_e10 * result)
  282. {
  283.   const double ay  = fabs(y);
  284.  
  285.   if(y == 0.0) {
  286.     result->val = 0.0;
  287.     result->err = fabs(dy * exp(x));
  288.     result->e10 = 0;
  289.     return GSL_SUCCESS;
  290.   }
  291.   else if(   ( x < 0.5*GSL_LOG_DBL_MAX   &&   x > 0.5*GSL_LOG_DBL_MIN)
  292.           && (ay < 0.8*GSL_SQRT_DBL_MAX  &&  ay > 1.2*GSL_SQRT_DBL_MIN)
  293.     ) {
  294.     const double ex = exp(x);
  295.     result->val  = y * ex;
  296.     result->err  = ex * (fabs(dy) + fabs(y*dx));
  297.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  298.     result->e10 = 0;
  299.     return GSL_SUCCESS;
  300.   }
  301.   else {
  302.     const double ly  = log(ay);
  303.     const double l10_val = (x + ly)/M_LN10;
  304.  
  305.     if(l10_val > INT_MAX-1) {
  306.       OVERFLOW_ERROR_E10(result);
  307.     }
  308.     else if(l10_val < INT_MIN+1) {
  309.       UNDERFLOW_ERROR_E10(result);
  310.     }
  311.     else {
  312.       const double sy  = GSL_SIGN(y);
  313.       const int    N   = (int) floor(l10_val);
  314.       const double arg_val = (l10_val - N) * M_LN10;
  315.       const double arg_err = dy/fabs(y) + dx + 2.0*GSL_DBL_EPSILON*fabs(arg_val);
  316.  
  317.       result->val  = sy * exp(arg_val);
  318.       result->err  = arg_err * fabs(result->val);
  319.       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  320.       result->e10 = N;
  321.  
  322.       return GSL_SUCCESS;
  323.     }
  324.   }
  325. }
  326.  
  327.  
  328. int gsl_sf_expm1_e(const double x, gsl_sf_result * result)
  329. {
  330.   const double cut = 0.002;
  331.  
  332.   if(x < GSL_LOG_DBL_MIN) {
  333.     result->val = -1.0;
  334.     result->err = GSL_DBL_EPSILON;
  335.     return GSL_SUCCESS;
  336.   }
  337.   else if(x < -cut) {
  338.     result->val = exp(x) - 1.0;
  339.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  340.     return GSL_SUCCESS;
  341.   }
  342.   else if(x < cut) {
  343.     result->val = x * (1.0 + 0.5*x*(1.0 + x/3.0*(1.0 + 0.25*x*(1.0 + 0.2*x))));
  344.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  345.     return GSL_SUCCESS;
  346.   } 
  347.   else if(x < GSL_LOG_DBL_MAX) {
  348.     result->val = exp(x) - 1.0;
  349.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  350.     return GSL_SUCCESS;
  351.   }
  352.   else {
  353.     OVERFLOW_ERROR(result);
  354.   }
  355. }
  356.  
  357.  
  358. int gsl_sf_exprel_e(const double x, gsl_sf_result * result)
  359. {
  360.   const double cut = 0.002;
  361.  
  362.   if(x < GSL_LOG_DBL_MIN) {
  363.     result->val = -1.0/x;
  364.     result->err = GSL_DBL_EPSILON * fabs(result->val);
  365.     return GSL_SUCCESS;
  366.   }
  367.   else if(x < -cut) {
  368.     result->val = (exp(x) - 1.0)/x;
  369.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  370.     return GSL_SUCCESS;
  371.   }
  372.   else if(x < cut) {
  373.     result->val = (1.0 + 0.5*x*(1.0 + x/3.0*(1.0 + 0.25*x*(1.0 + 0.2*x))));
  374.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  375.     return GSL_SUCCESS;
  376.   } 
  377.   else if(x < GSL_LOG_DBL_MAX) {
  378.     result->val = (exp(x) - 1.0)/x;
  379.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  380.     return GSL_SUCCESS;
  381.   }
  382.   else {
  383.     OVERFLOW_ERROR(result);
  384.   }
  385. }
  386.  
  387.  
  388. int gsl_sf_exprel_2_e(double x, gsl_sf_result * result)
  389. {
  390.   const double cut = 0.002;
  391.  
  392.   if(x < GSL_LOG_DBL_MIN) {
  393.     result->val = -2.0/x*(1.0 + 1.0/x);
  394.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  395.     return GSL_SUCCESS;
  396.   }
  397.   else if(x < -cut) {
  398.     result->val = 2.0*(exp(x) - 1.0 - x)/(x*x);
  399.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  400.     return GSL_SUCCESS;
  401.   }
  402.   else if(x < cut) {
  403.     result->val = (1.0 + 1.0/3.0*x*(1.0 + 0.25*x*(1.0 + 0.2*x*(1.0 + 1.0/6.0*x))));
  404.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  405.     return GSL_SUCCESS;
  406.   } 
  407.   else if(x < GSL_LOG_DBL_MAX) {
  408.     result->val = 2.0*(exp(x) - 1.0 - x)/(x*x);
  409.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  410.     return GSL_SUCCESS;
  411.   }
  412.   else {
  413.     OVERFLOW_ERROR(result);
  414.   }
  415. }
  416.  
  417.  
  418. int
  419. gsl_sf_exprel_n_e(const int N, const double x, gsl_sf_result * result)
  420. {
  421.   if(N < 0) {
  422.     DOMAIN_ERROR(result);
  423.   }
  424.   else if(x == 0.0) {
  425.     result->val = 1.0;
  426.     result->err = 0.0;
  427.     return GSL_SUCCESS;
  428.   }
  429.   else if(fabs(x) < GSL_ROOT3_DBL_EPSILON * N) {
  430.     result->val = 1.0 + x/(N+1) * (1.0 + x/(N+2));
  431.     result->err = 2.0 * GSL_DBL_EPSILON;
  432.     return GSL_SUCCESS;
  433.   }
  434.   else if(N == 0) {
  435.     return gsl_sf_exp_e(x, result);
  436.   }
  437.   else if(N == 1) {
  438.     return gsl_sf_exprel_e(x, result);
  439.   }
  440.   else if(N == 2) {
  441.     return gsl_sf_exprel_2_e(x, result);
  442.   }
  443.   else {
  444.     if(x > N && (-x + N*(1.0 + log(x/N)) < GSL_LOG_DBL_EPSILON)) {
  445.       /* x is much larger than n.
  446.        * Ignore polynomial part, so
  447.        * exprel_N(x) ~= e^x N!/x^N
  448.        */
  449.       gsl_sf_result lnf_N;
  450.       double lnr_val;
  451.       double lnr_err;
  452.       double lnterm;
  453.       gsl_sf_lnfact_e(N, &lnf_N);
  454.       lnterm = N*log(x);
  455.       lnr_val  = x + lnf_N.val - lnterm;
  456.       lnr_err  = GSL_DBL_EPSILON * (fabs(x) + fabs(lnf_N.val) + fabs(lnterm));
  457.       lnr_err += lnf_N.err;
  458.       return gsl_sf_exp_err_e(lnr_val, lnr_err, result);
  459.     }
  460.     else if(x > N) {
  461.       /* Write the identity
  462.        *   exprel_n(x) = e^x n! / x^n (1 - Gamma[n,x]/Gamma[n])
  463.        * then use the asymptotic expansion
  464.        * Gamma[n,x] ~ x^(n-1) e^(-x) (1 + (n-1)/x + (n-1)(n-2)/x^2 + ...)
  465.        */
  466.       double ln_x = log(x);
  467.       gsl_sf_result lnf_N;
  468.       double lg_N;
  469.       double lnpre_val;
  470.       double lnpre_err;
  471.       gsl_sf_lnfact_e(N, &lnf_N);    /* log(N!)       */
  472.       lg_N  = lnf_N.val - log(N);       /* log(Gamma(N)) */
  473.       lnpre_val  = x + lnf_N.val - N*ln_x;
  474.       lnpre_err  = GSL_DBL_EPSILON * (fabs(x) + fabs(lnf_N.val) + fabs(N*ln_x));
  475.       lnpre_err += lnf_N.err;
  476.       if(lnpre_val < GSL_LOG_DBL_MAX - 5.0) {
  477.         int stat_eG;
  478.     gsl_sf_result bigG_ratio;
  479.     gsl_sf_result pre;
  480.     int stat_ex = gsl_sf_exp_err_e(lnpre_val, lnpre_err, &pre);
  481.         double ln_bigG_ratio_pre = -x + (N-1)*ln_x - lg_N;
  482.     double bigGsum = 1.0;
  483.     double term = 1.0;
  484.     int k;
  485.     for(k=1; k<N; k++) {
  486.       term *= (N-k)/x;
  487.       bigGsum += term;
  488.     }
  489.     stat_eG = gsl_sf_exp_mult_e(ln_bigG_ratio_pre, bigGsum, &bigG_ratio);
  490.     if(stat_eG == GSL_SUCCESS) {
  491.           result->val  = pre.val * (1.0 - bigG_ratio.val);
  492.       result->err  = pre.val * (2.0*GSL_DBL_EPSILON + bigG_ratio.err);
  493.       result->err += pre.err * fabs(1.0 - bigG_ratio.val);
  494.       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  495.       return stat_ex;
  496.     }
  497.     else {
  498.       result->val = 0.0;
  499.       result->err = 0.0;
  500.       return stat_eG;
  501.     }
  502.       }
  503.       else {
  504.         OVERFLOW_ERROR(result);
  505.       }
  506.     }
  507.     else if(x > -10.0*N) {
  508.       return exprel_n_CF(N, x, result);
  509.     }
  510.     else {
  511.       /* x -> -Inf asymptotic:
  512.        * exprel_n(x) ~ e^x n!/x^n - n/x (1 + (n-1)/x + (n-1)(n-2)/x + ...)
  513.        *             ~ - n/x (1 + (n-1)/x + (n-1)(n-2)/x + ...)
  514.        */
  515.       double sum  = 1.0;
  516.       double term = 1.0;
  517.       int k;
  518.       for(k=1; k<N; k++) {
  519.         term *= (N-k)/x;
  520.     sum  += term;
  521.       }
  522.       result->val = -N/x * sum;
  523.       result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  524.       return GSL_SUCCESS;
  525.     }
  526.   }
  527. }
  528.  
  529.  
  530. int
  531. gsl_sf_exp_err_e(const double x, const double dx, gsl_sf_result * result)
  532. {
  533.   const double adx = fabs(dx);
  534.  
  535.   /* CHECK_POINTER(result) */
  536.  
  537.   if(x + adx > GSL_LOG_DBL_MAX) {
  538.     OVERFLOW_ERROR(result);
  539.   }
  540.   else if(x - adx < GSL_LOG_DBL_MIN) {
  541.     UNDERFLOW_ERROR(result);
  542.   }
  543.   else {
  544.     const double ex  = exp(x);
  545.     const double edx = exp(adx);
  546.     result->val  = ex;
  547.     result->err  = ex * GSL_MAX_DBL(GSL_DBL_EPSILON, edx - 1.0/edx);
  548.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  549.     return GSL_SUCCESS;
  550.   }
  551. }
  552.  
  553.  
  554. int
  555. gsl_sf_exp_err_e10_e(const double x, const double dx, gsl_sf_result_e10 * result)
  556. {
  557.   const double adx = fabs(dx);
  558.  
  559.   /* CHECK_POINTER(result) */
  560.  
  561.   if(x + adx > INT_MAX - 1) {
  562.     OVERFLOW_ERROR_E10(result);
  563.   }
  564.   else if(x - adx < INT_MIN + 1) {
  565.     UNDERFLOW_ERROR_E10(result);
  566.   }
  567.   else {
  568.     const int    N  = (int)floor(x/M_LN10);
  569.     const double ex = exp(x-N*M_LN10);
  570.     result->val = ex;
  571.     result->err = ex * (2.0 * GSL_DBL_EPSILON * (fabs(x) + 1.0) + adx);
  572.     result->e10 = N;
  573.     return GSL_SUCCESS;
  574.   }
  575. }
  576.  
  577.  
  578. /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
  579.  
  580. #include "eval.h"
  581.  
  582. double gsl_sf_exp(const double x)
  583. {
  584.   EVAL_RESULT(gsl_sf_exp_e(x, &result));
  585. }
  586.  
  587. double gsl_sf_exp_mult(const double x, const double y)
  588. {
  589.   EVAL_RESULT(gsl_sf_exp_mult_e(x, y, &result));
  590. }
  591.  
  592. double gsl_sf_expm1(const double x)
  593. {
  594.   EVAL_RESULT(gsl_sf_expm1_e(x, &result));
  595. }
  596.  
  597. double gsl_sf_exprel(const double x)
  598. {
  599.   EVAL_RESULT(gsl_sf_exprel_e(x, &result));
  600. }
  601.  
  602. double gsl_sf_exprel_2(const double x)
  603. {
  604.   EVAL_RESULT(gsl_sf_exprel_2_e(x, &result));
  605. }
  606.  
  607. double gsl_sf_exprel_n(const int n, const double x)
  608. {
  609.   EVAL_RESULT(gsl_sf_exprel_n_e(n, x, &result));
  610. }
  611.